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ABSTRACT 

Differencing  operators  of  arbitrarily  high  order  can  be  constructed  by  in¬ 
terpolating  a  polynomial  through  a  set  of  data  followed  by  differentiation  of 
this  polynomial  and  finally  evaluation  of  the  polynomial  at  the  point  where  a 
derivative  approximation  is  desired.  Furthermore,  the  interpolating  polyno¬ 
mial  can  be  constructed  from  algebraic,  trigonometric,  or,  perhaps  exponen¬ 
tial  polynomials.  This  paper  begins  with  a  comparison  of  such  differencing 
operator  construction.  Next,  the  issue  of  proper  grids  for  high  order  poly¬ 
nomials  is  addressed.  Finally,  an  adaptive  numerical  method  is  introduced 
which  adapts  the  numerical  grid  and  the  order  of  the  differencing  operator 
depending  on  the  data.  The  numerical  grid  adaptation  is  performed  on  a 
Chebyshev  grid.  That  is,  at  each  level  of  refinement  the  grid  is  a  Chebyshev 
grid  and  this  grid  is  refined  locally  based  on  wavelet  analysis. 
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1  Introduction 


One  can  argue  that  high  order  numerical  methods  are  appropriate  for  prob¬ 
lems  in  1)  the  direct  numerical  simulation  of  turbulence,  see  [17],  2)  flows 
with  shocks  and  nonlinear  physics,  see  [8]  and  3)  flows  with  smooth  propagat¬ 
ing  structures  such  as  those  encountered  in  aeroacoustics.  Assertion  number 
3)  is  based  on  convergence  properties  of  the  hp-refinement  method  in  finite 
elements,  see  [20],  [1],  [11],  in  which  convergence  is  very  fast  for  high  order 
polynomials  as  long  as  the  function  at  hand  is  smooth.  In  addition,  high 
order  methods  are  more  efficient  for  long  time  integration  of  unsteady  flow 
problems,  see  [18]. 

This  paper  introduces  a  numerical  method  which  combines  very  high 
order  differencing  with  a  wavelet-based  grid  and  order  selection  mechanism. 
Here  very  high  order  differencing  will  be  schemes  of  order  greater  than  or 
equal  to  8,  i.e.,  perhaps  16,  20,  or  maybe  order  32.  Such  high  orders  of 
accuracy  can  produce  solutions  which  are  very  close  to  those  produced  by 
spectral  methods.  See  [3]  for  a  spectral  method  on  arbitrary  grids. 

The  numerical  method  introduced  here  is  named  the  Wavelet- Optimized 
Finite  Difference  method  2,  or  WOFD2.  WOFD2  is  an  extension  of  WOFD, 
see  [15],  [9],  in  which  wavelets  are  used  for  grid  refinement,  see  [16],  for  finite 
difference  schemes.  The  new  method  extends  this  scheme  to  very  high  orders 
and  also  adapts  the  order  of  accuracy  depending  on  the  data. 

Let  us  begin  by  studying  various  manners  in  which  high  order  difference 
operators  can  be  constructed. 


2  Generating  Difference  Equations 

Given  a  vector  of  N  numbers  /  how  can  we  get  an  approximate  value  of  the 
derivative  f  at  the  i  —  th  point  and  how  good  will  this  approximate  value  be. 
Generally  speaking,  the  more  elements  around  the  i  —  th  point  of  /  that  are 
used  to  approximate  /'  the  better  the  approximation  will  be.  Common  finite 
difference  formulas  are  found  by  fitting  a  algebraic  polynomial  of  degree  q  lo¬ 
cally  around  the  i—th  point  of  a  vector  /  of  evenly-spaced  elements  to  obtain 
difference  approximations  of  accuracy  q  — I .  This  section  will  generalize  this 
concept  to  find  the  difference  equations  of  arbitrary  accuracy  on  arbitrary 
grids  using  algebraic,  trigonometric,  cosine  and  exponential  polynomials.  As 
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special  cases,  one  can  obtain  all  the  usual  finite  difference  formulas  as  well  as 
the  Fourier  collocation  and  Cbebyshev  collocation  spectral  differential  ma¬ 
trices. 

Two  methods  of  generating  the  differencing  coefficients  will  be  introduced. 
The  first  method  explains  how  to  set  up  a  system  of  equations  which  will  have 
as  a  solution  the  differencing  coefficients.  The  second  method  is  the  deriva¬ 
tion  of  differencing  coefficients  by  interpolation.  It  is  this  second  method 
which  is  used  throughout  the  paper  for  the  actual  generation  of  difference 
equations. 


2.1  Setting  up  a  Linear  System 

The  problem  is  to  find  a  set  of  coefficients  {fjt}  which  combines  the  raw  data 
in  a  vector  /  to  provide  an  approximation  to  a  derivative: 

k=righi 

Z)  ^kfixk).  (1) 

k-left 

If  we  require  the  above  equation  to  be  exact  for  polynomials,  algebraic,trigonometric, 
cosine  or  exponential,  then  a  linear  system  of  equations  can  be  solved  to  find 
an  appropriate  set  of  differencing  coefficients  {r^}.  Let  b{x)  denote  a  funda¬ 
mental  basis  element  from  which  a  basis  can  generated  by  taking  powers  of 
6(x):  b{x)  =  X,  b{x)  =  e‘®,  b{x)  =  cos(x),  or  b(x)  =  e®.  That  is,  we  require 
that  the  derivative  be  exact  up  to  a  given  order  N  on  the  numerical  grid. 

The  system  of  equations  to  be  solved  for  a  centered  differentiation  stencil  is 
as  follows: 

=  E  n{b{xi^,)r,  (2) 

k=-L 

iFrom  this  equation  one  can  generate  a  system  of  equations  with  N  re¬ 
quirements,  that  N  functions  be  differentiated  exactly,  and  N  degrees  of 
freedom,  the  N  differencing  coefficients  rjt.  If  one  is  near  a  boundary,  then 
the  stenciled  is  biased.  Since  this  type  of  system  is  well-known  for  algebraic 
polynomials,  an  example  for  the  less  well-known  trigonometric  polynomials 
will  be  given. 

Consider  a  trigonometric  polynomial  on  a  3  point  centered  stencil.  The 
first  equation  simply  requires  that  the  derivative  of  a  constant  be  zero: 

0  =  r_i  -I-  To  +  ri.  (3) 
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Note  that  is  the  same  equation  as  for  algebraic  polynomials  since  (a:)°  = 
(e’®)°.  The  next  two  equations  come  from  requiring  that  the  n  =  1  mode  is 
differentiated  exactly: 


ie’*®  =  r_i  +  roe’®°  +  ri  .  (4) 

One  now  obtains  the  two  equations  from  equating  the  real  and  imaginary 
parts.  These  three  equations  can  be  solved  for  the  three  coefficients  r_i,  ro,  ri. 

Similarly,  one  can  find  the  coefficients  for  higher  order  schemes  by  requir¬ 
ing  that  more  modes  be  differentiated  exactly.  Note  that  no  restrictions  were 
placed  on  the  grid.  Differencing  formulas  can  be  found  on  arbitrary  grids  as 
easily  as  they  can  be  found  on  uniform  grids.  Also,  note  that  the  Fourier 
spectral  differentiation  matrix  can  be  found  from  the  above  procedure  by 
requiring  that  the  grid  be  uniform  and  that  the  differencing  formulas  have 
maximum  accuracy  on  a  given  grid.  That  is,  if  one  is  working  on  a  grid 
of  size  33  then  require  that  the  first  16  modes  and  the  zero-th  mode  are 
differentiated  exactly. 

2.2  Interpolation 

A  second  approach,  and  the  one  used  in  this  paper,  is  to  generate  differencing 
coefficients  by  first  interpolating  a  polynomial  through  a  set  of  data,  followed 
by  differentiation  of  this  polynomial  and  evaluated  at  a  grid  point. 

The  main  reason  that  differentiation  was  studied  with  a  variety  of  types  of 
differentiation  operators  was  to  find  out  if  there  was  any  advantage  to  using, 
say,  trigonometric  polynomials  to  differentiate  as  opposed  to  algebraic  poly¬ 
nomials  when  the  function  to  be  differentiated  was  for  example  a  Gaussian 
pulse.  It  seemed  like  an  appropriate  study  to  undertake  given  the  current 
research  activity  in  the  area  or  aeroacoustics  where  one  is  often  confronted 
with  the  need  to  computationally  propagate  some  type  of  wave  motion.  The 
thought  was  that  perhaps  trigonometric  polynomials  might  have  some  ad¬ 
vantage  at  propagating  wave  motion  over  the  more  common  algebraic  poly¬ 
nomials.  One  of  the  conclusions  of  this  section  is  that  there  is  no  advantage 
and  that  one  should  simply  use  algebraic  polynomials  for  the  generation  of 
differencing  equations.  In  fact,  the  only  important  issues  involved  with  ob¬ 
taining  approximate  derivatives  is  the  order  of  the  finite  difference  operator 
and  the  density  of  the  numerical  grid. 
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The  most  important  reference  for  this  section  is  [7].  The  following  four 
subsections  will  cite  the  interpolation  formulas  for  the  four  types  of  interpo¬ 
lation,  and  hence  diiferentiation,  considered  in  this  section. 

2.2.1  Algebraic  Polynomials 

Interpolation  with  algebraic  polynomials  is  probably  the  most  common  form 
of  interpolation,  and  it  is  from  this  type  of  interpolation  that  common  uni¬ 
form  grid  finite  difference  methods  can  be  formd.  Using  the  following  for¬ 
mula  one  can  find  the  finite  difference  coefficients  for  an  arbitrary  grid  and 
of  arbitrary  order.  One  simply  fits  the  polynomial  to  the  data,  followed  by 
differentiation  of  the  polynomial,  and  finally  one  evaluates  the  polynomial 
at  the  point  of  interest.  The  well-known  Lagrange  interpolation  formula  for 
algebraic  interpolation  is, 

=  n  n  (5) 

k=:Ojkj^j  k=Oyk^j 

Aj{xk)  =  Sjk  For  given  values  too,  the  polynomial 

Pn{x)  =  E  fJ^kAk{x).  (6) 

fc=0 

in  Pn  and  takes  on  these  values  at  the  points  a;,-: 

Pn{xk)  =  m,  (7) 


for  A:  =  0, 1, ...,  n. 

2.2.2  Trigonometric  Polynomials 

As  seen  from  the  previous  section,  one  can  also  generate  difference  operators 
by  using  trigonometric  functions  as  the  fundamental  interpolation  elements. 
The  following  is  the  appropriate  Lagrange-type  interpolation  formula,  see  [7]; 
For  —7r<xo<xi<  ...  <  X2n  <  then 

2n  ^  2n  2 

Tj{x:)=  n  sin-(a:-a;i;)/  JJ  sm-{xj-Xk).  (8) 

k=0,k^j  k=0,k^j 
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The  function, 
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T{x)  =  Y,WkTk{x)  (9) 

Jfe=0 

is  the  unique  solution  of  the  interpolation  problem, 

T(xfc)  =  tiJfc,  (10) 

for  k  =  0, 1,..., 2n.  Again,  one  can  derive  finite  difference  coefficients  by 
interpolating  to  a  function,  followed  by  differentiation  of  the  interpolation 
polynomial  and  evaluation  at  the  point  of  interest.  The  following  section  will 
prove  that  such  difference  equations  obey  order  properties  just  as  the  usual 
difference  equations  derived  from  algebraic  polynomials  do. 

2.2.3  Cosine  Polynomials 

The  comparable  Lagrange-style  interpolation  formula  for  cosine  polynomials 
is  the  following,  see  [7]. 

Given  n  +  1  distinct  points  0  <  xo  <  xi  <  ...  <  x„  <  -k.  Set 

n  n 

n  (cosx  —  COSXfc)/  JJ  (coSXj  —  COSXfc).  (11) 

Then  Cj  is  a  cosine  polynomial  of  order  <  n,  Cj{x)  =  jyk=o^kcos{kx),  for 
which  Cj{xk)  =  Sjk-  Given  n+1  distinct  values  wq,  wi, ...,  Wn  there  is  a  unique 
cosine  polynomial  of  order  <  n,  C(x),  for  which  C{xk)  =  Wk,  fc  =  0, 1, ...,  n. 

n 

H  WkCk{x).  (12) 

Note  that  £(7(x)|o,x  =  0  since  ^C'fc(x)|o,^  =  0  for  all  k.  For  this  rea¬ 
son,  difference  operators  based  on  cosine  polynomials  will  not  be  considered 
in  general,  but  will  be  compared  in  a  later  section  to  Chebyshev  spectral 
methods.  As  above,  the  differencing  coefficients  are  found  by  first  fitting 
the  trigonometric  polynomial  to  the  data,  followed  by  differentiation  of  the 
polynomial  and  finally  evaluation  at  the  point  of  interest. 
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2.2.4  Exponential  Polynomials 

The  final  polynomial  to  be  tested  is  the  exponential  polynomial, 

n  («'-«**)/  n  (e*'-e'‘), 

where  the  interpolation  polynomial  is, 

n 

=  £  '>^kEk{x), 
k=0 

and 

E{xk)  =  Wk. 


(13) 


(14) 

(15) 


2.3  Truncation  Error  and  Differentiation  Accuracy 

The  purpose  of  this  section  is  to  illustrate  algebraically  that  one  obtains  dif¬ 
ferentiation  order-of- accuracy  properties  for  all  four  types  of  differentiation 
operators  which  are  similar  to  the  standard  order-of-accuracy  properties  ob¬ 
tained  with  the  usual  algebraic  interpolation.  In  short,  if  one  interpolates 
an  N  order  polynomial  then  one  obtains  a  reduction  in  differentiation  error 
of  (|)^  when  the  density  of  the  grid  is  doubled.  This  order  of  accuracy  is 
obtained  regardless  of  the  type  of  polynomials  which  are  used. 

Recall  that  the  remainder  for  algebraic  polynomial  interpolation  is,  see 

p]. 

/(,)  _  --“)/("-■>«),  (w) 

where  ^  lies  between  the  smallest  and  the  largest  x,-.  The  following  sec¬ 
tion  will  show  that  a  similar  expression  can  be  obtained  for  any  polynomial 
constructed  from  powers  of  a  given  function. 

2.3.1  Truncation  Error  for  Interpolation  by  Powers 

There  are  some  subtle  issues  concerning  a  general  proof  of  truncation  error 
and  accuracy  for  interpolation  by  a  polynomial  constructed  from  the  powers 
of  a  general  function  g{x),  see  [4].  The  following  demonstration  will  illustrate 
the  essentual  algebraic  steps  that  one  follows  to  obtain  accurcy  while  avoiding 
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the  subtle  issues.  In  short,  let  the  polynomial  element  g{x)  and  the  function 
to  be  approximated  /(x)  be  “well-behaved”. 

n 

A;=0 


be  the  polynomial  which  interpolates  f{x)  at  xo,a;i,...,a:„,  p{xi) 
then, 

f(x)-p(x)  = 


where 


^(x)  =  (^(x)  -  g{xo)){g{x)  -  g(xi))...{g{x)  -  g{xn)), 
and  where  ^  lies  between  the  smallest  and  the  largest  x^. 


(17) 

(18) 


Demonstration  of  Truncation  Error:  Note  that  much  of  this  demon¬ 
stration  is  the  same  as  that  which  can  be  found  in  a  standard  numerical 
analysis  text  for  the  remainder  term  in  algebraic  interpolation,  see  [5]. 
Define  H{z)  such  that 


H{z)  =  f{z)  -  p{z)  -  R{x)(f>{z),  (19) 


where  R{x)  is  defined  such  that  H{x)  =  0.  Note  that  H{xi)  =  0,  for  i  = 
0,  ...,n  since  p(x,)  =  f{xi)  and  ^(x,)  =  0.  i,From  Rolle’s  theorem  it  follows 
that  there  exists  a  point  ^  in  the  interval  defined  by  the  smallest  and  largest 
Xj’s  such  that  =  0.  This  implies, 


R{x) 


/(-»)({)  _p(.w)(^) 


(20) 


Now  put  this  back  into  the  expression  for  H{z)  and  set  .2  =  x  to  get. 


/(x)  -p(x) 


/(n+l)(^)-yn-H)(^) 

^(«+l)(^) 


4>{x). 


(21) 


This  is  the  desired  expression.  / /. 
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Note  that  in  the  above  demonstration  that  if  the  polynomial  is  algebraic 
that  =  0  and  =  (n  +  1)!,  but  for  a  general  g{x)  these 

two  functions  are  just  a  measure  of  the  smoothness  of  the  basic  interpolation 
element,  g(x).  still  remains  a  measure  of  the  smoothness  of  the 

function  one  is  interpolating  to,  and  (f>{x)  is  a  function  dependent  on  the  grid 
distribution. 

2.3.2  Differentiation  Accuracy 

The  primary  interest  here  is  to  understand  the  behavior  of  the  derivative 
operators  derived  from  the  various  types  of  interpolation  outlined  above  as 
the  grid  is  refined.  That  is, 

f{x)-p'{x)  =  Q{On^).  (22) 

where  Q{^)  =  ,  and  (f>'{x)  will  dictate  the  behavior  as  the 

grid  is  refined.  It  will  be  shown  that  the  behavior  of  is  essentially 
independent  of  the  basic  interpolation  element  g{x)  and  depends  only  on  the 
order  of  the  interpolation. 

Demonstration  of  Accuracy: 

Let  h  denote  the  smallest  spacing  in  the  numerical  grid,  then 

<f>'{xo)  =  Chr  +  h.o.t.  (23) 

where  n  is  the  highest  power  in  the  interpolation  polynomial,  and  the  point 
xq  is  an  arbitrary  grid  point  inside  the  interpolation  stencil. 

Demonstration:  First  of  all, 

(f>\xo)  =  g'{xo){g{xo)  -  g{x-i){g{xo)  -  g{x2))...{g{xa:)  -  g{xn))-  (24) 

Expand  about  zero  the  function  g{x), 

g{x)  =  if(0)  +  g'{^)x  +  5r"(0)x^/2  +  ...  (25) 

and  examine  the  difference  g{xo)  —  g{x\) 

g{xo)-g{xx)  ^  g'mxo-x^)^g"{Q)l2{xl-xl)+g"\Q)l^{xl-x\)  +  ...  (26) 
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Without  loss  of  generality  let  one  of  the  points  be  zero,  say  xq  =  0,  to  get, 
s(0)  -  =  s'(0)(-x0  +  s"(0)/2(-x;)  +  ff"'(0)/6(-x?)  +  ....  (27) 


or, 

j(0)-s(x.)  =  -x.(i;^xr'),  (28) 

and  one  can  see  that  the  first  term  in  the  difference  is  linear.  If  xo  is  not  zero 
then  one  obtains  the  factorization, 

X®  -  j/5  =  (x  -  (29) 

»=o 


and  hence, 


oo  771-1 

di^o)  -  g{xi)  =  (xo  -  Xi)(f^  ~7|  (30) 

It  should  be  clear  that  the  first  term  in  the  difference  g(xo)  —  g(xi)  is  the 
linear  term  and  that  doubling  the  grid  such  that  between  every  two  points 
another  point  is  placed  is,  therefore,  halves  the  distance  to  a  first  order 
approximation.  For  each  of  the  differences  Xj  —  xj  there  exists  a  constant  c,j 
such  that  the  difference, 

X,'  “■  Xj  —  (31 ) 

can  be  expressed  in  terms  of  the  smallest  grid  spacing  h.  Let  C  =  11^=1  coj 
then  it  is,  therefore,  clear  that 

<f>'{xo)  =  Ch^  +  h.o.t.  (32) 


//. 

Consider  the  following  special  cases  which  includes  all  the  polynomial 
types  discussed  above: 

•  ^f(x)  =  X,  5r^’"^(0)  =  0,  for  m  1 

•  g{x)  =  e®,  5r(™)(0)  =  1,  Vm 
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For  a  simple  illustration,  consider  the  following  two  examples  of  algebraic 
and  exponential  interpolation. 

Algebraic  Interpolation 

Consider  the  simple  case  of  interpolating  an  algebraic  quadratic  polyno¬ 
mial  P2{x)  to  a  function  f{x)  at  the  grid  points  rco  <  a;i  <  x^:  P2{xi)  =  f{xi), 
i  =  0, 1, 2.  The  remainder  term  for  some  Xo<(  <  X2,  is 

f{x)  -  p2{x)  =  {x-  xo)(a;  -  xi){x  -  (33) 

Now,  differentiate  and  evaluate  at  x  =  xi  to  get, 

fix,)  -  =  (xi  -  (34) 

where  h  =  Xi  —  xq  =  X2  —  Xi.  If  the  grid  is  evenly-spaced  then  the  differences 
(xi  —  Xj)  are  some  integer  multiple  of  the  smallest  difference  which  one  can 
denote  by  h.  If  one  doubles  the  number  of  grid  points  then  each  of  the 
distances  (xj  —  Xj)  becomes  half  as  large  and  the  accuracy  for  this  quadratic 
example  will  be  2. 

The  General  Statement  for  Algebraic  Interpolation 

In  general,  one  can  expect  that  algebraic  polynomial  interpolation  with 
a  polynomial  of  order  n,  Pn(x),  will  produce  a  differencing  operator  of  order 
also  of  order  n.  This  can  be  seen  from  the  portion  of  the  truncation  error 
which  depends  on  the  grid  distribution: 

n(a;  -  Xi).  (35) 

i=0 

This  product  contains  n  -|- 1  terms.  After  differentiation  and  evaluation  at  a 
point  Xk  the  product  will  contain  n  terms, 

Jl(xfc-Xi).  (36) 

i=0 
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When  the  grid  density  is  doubled  by  adding  a  point  between  every  two  points, 
then  each  distance  {xk  —  x,)  will  decrease  by  a  factor  of  2.  The  product  of 
these  factors  will  be  a  multiple  of  (I)"  and  hence  the  accuracy  of  n  is  achieved. 

Exponential  Interpolation 

If,  on  the  other  hand,  P2(x)  is  now  an  exponential  polynomial  then  after 
differentiation  and  evaluation  one  gets, 

/'(x.)  -  pi(x.)  =  e«(e-  -  e“)(e«  -  (37) 

Without  loss  of  generality,  let  xx  =  0.  Then  the  product, 

/'(O)  -  ?2(0)  =  (38) 

becomes 

/'(O)  —^2(0)  ~  (^0  +  a:o/2  +  ..•)(^2  +  +  ••O?  (39) 

and  one  can  see  that  if  the  distance  to  zero  is  halved  and  hence  xq  and  X2  are 
divided  by  2  that  the  leading  order  terms  in  each  of  the  above  parenthesis 
dictates  that  the  error  will  be  reduced  by  4  to  a  first  order  approximation. 

Therefore,  one  can  expect  the  accuracy  to  behave  as  it  does  for  algebraic 
polynomials.  That  is,  doubling  the  grid  points  will  decrease  the  error  by 
(1  /2)^  to  first  order. 

The  General  Statement  for  Exponential  Interpolation 

In  general,  one  can  expect  that  exponential  polynomial  interpolation  with 
a  polynomial  of  order  n,  Pn{x),  wiU  produce  a  differencing  operator  of  order 
also  of  order  n  which  is  the  same  result  as  for  algebraic  interpolation.  As  can 
be  seen  from  the  above  example,  differentiation  and  evaluation  at  a  point  Xk 
will  produce  a  leading  order  term  which  is  a  product  of  n  terms, 

n-1 

n(a:fc-a:i),  (40) 

z=0 

and  accuracy  of  order  n  is  achieved. 
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Trig 

Err 

Err 

Ratio 

Exp 

Err 

Err 

Ratio 

mB 

4.10  *  10-^ 

^1 

64 

EH 

mm 

29.3 

128 

9.54  *  10-^2 

aMliH 

29.8 

9.04  =*=  10~^® 

EH 

512 

9.70  *  10"^*^ 

210.0 

8.95  *  10-^« 

210.0 

2.17  *  10“^^ 

210.0 

Table  1:  A  Comparison  of  Errors  for  Various  Types  of  Differentiation  Oper¬ 
ators 


2.4  A  Numerical  Check  of  Accuracy 

This  subsection  will  verify  that  the  differentiation  operators  which  are  gener¬ 
ated  from  algebraic,  trigonometric,  and  exponential  polynomials  all  exhibit 
the  same  order  property,  which  depends  only  on  the  order  of  the  polynomial 
interpolation.  The  difference  operators  will  be  tested  on  the  function, 

^w=2Ti(S) 

defined  on  [0,  tt].  f{x)  is  chosen  because  it  is  periodic  but  not  exactly  a 
trigonometric  or  algebraic  polynomial.  Table  (1)  illustrates  the  order  prop¬ 
erty  All  of  the  errors  are  L2. 


A  general  question  arises,  is  there  any  advantage  to  using,  say,  a  trigono¬ 
metric  polynomial  for  the  generation  of  difference  equations  over,  say,  the 
usual  algebraic  polynomial?  From  this  study  the  answer  appears  to  be  no. 
The  essence  depends  on  the  ability  of  the  interpolating  polynomial  to  locally 
approximate  the  function  at  hand.  If  the  function  at  hand  is  not  exactly  a 
trigonometric  or  algebraic  polynomial,  as  is  likely,  then  there  is  no  advantage 
for  either  approach.  Such  an  issue  is  important  when  one  is  considering  the 
propagation  of,  say,  a  pulse  in  the  application  of  aeroacoustics.  A  pulse  will 
locally  be  neither  a  algebraic  or  trigonometric  polynomial.  In  short,  a  wave. 
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or  pulse,  can  not  be  propagated  accurately  if  it  can  not  first  be  differentiated 
accurately,  and  it  can  not  be  differentiated  accurately  if  it  can  not  first  be 
approximated  accurately. 

3  High  Order  Methods 

Spectral  collocation  methods  are  often  given  the  probable  misnomer  of  “in¬ 
finitely  accurate”.  In  a  manner  consistent  with  finite  difference  methods, 
the  accuracy  of  spectral  collocation  methods  will  be  assigned  the  accuracy 
of  iV  —  1  when  applied  on  a  grid  of  N  points. 

This  section  will  begin  by  connecting  spectral  collocation  methods  to  fi¬ 
nite  difference  methods.  That  is,  spectral  methods  will  be  viewed  from  the 
point  of  view  of  the  maximum  finite  difference  method  on  a  given  grid.  Fol¬ 
lowing  the  comments  on  this  connection,  a  case  will  be  made  for  applying  very 
high  order  algebraically  generated  finite  difference  operators  on  Cheybshev 
grids  or,  equivalently,  applying  very  high  order  cosine  polynomial  generated 
finite  difference  operators  on  a  uniform  grid.  This  second  process  of  using 
cosine  polynomials  will  require  a  mapping  of  the  independent  variable  from, 
say,  X  to  cos(x),  but  is  exactly  equal  to  applying  algebraic  polynomials  on 
Chebyshev  grids. 

3.1  Spectral  Collocation  =  Maximum  Order  Finite 
Difference 

On  a  numerical  grid  of  N  points  one  can  fit  a  polynomial  with  N  degrees-of- 
freedom  through  all  of  the  data.  If  this  polynomial  is  algebraic  and  if  the  grid 
distribution  is  Chebyshev,  Xi  =  cos(^),  then  one  can  build  the  Chebyshev 
collocation  differentiation  matrix.  On  the  other  hand,  if  the  polynomial  is 
trigonometric  and  if  the  grid  is  uniformly  distributed  then  one  can  build 
the  Fourier  spectral  collocation  differentiation  matrix.  One  can,  therefore, 
define  in  the  physical  space  a  spectral  method  to  be  a  method  which  uses  the 
maximum  size  polynomial  for  approximation  and  differentiation  for  a  given 
grid  size. 

Another  way  to  see  this  is,  suppose  one  has  a  numerical  grid  of  16  points 
and  a  4-th  order  difference  operator  on  a  5  point  stencil.  Now  reduce  the 
number  of  grid  points  to  8.  The  difference  operator  is  still  4-th  order.  Now 
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Grid 

sin 

L2 

Pts 

Freq 

Error 

9 

1.0 

1.6710-28 

9 

1.5 

8.7110-^ 

9 

2.0 

5.6810-29 

9 

2.5 

1.6010" 

9 

3.0 

2.2410-28 

9 

3.5 

2.7910" 

9 

4.0 

8.010-28 

9 

4.5 

2.9110" 

9 

5.0 

3.0510" 

Table  2:  Fourier  Spectral  Collocation  Applied  to  sin’s 


reduce  the  number  of  grid  points  to  5.  The  difference  operator  is  still  4- 
th  order  accurate.  This  is  spectral  accuracy.  That  is,  spectral  accuracy  of 
collocation  methods  on  finite  grids  is  iV  —  1  where  N  is  the  number  of  grid 
points. 

3.1.1  A  Numerical  Check 

Let  us  consider  the  above  statements  numerically.  A  Fourier  collocation  spec¬ 
tral  method  on  a  grid  of  9  points  is  a  differencing  mechanism  with  exactly 
9  degrees-of-freedom,  and  hence,  is  designed  to  differentiate  exactly  9  func¬ 
tions  exactly:  1,  cos{kx),  and  sin(A;a:),  for  k  —  1,2, 3,4.  Table  (2)  is  meant 
to  illustrate  two  points:  i)  the  maximum  frequency  which  is  differentiated 
exactly  is  4.0,  and  ii)  the  poor  performance  on  non-integer  frequencies. 

Likewise,  a  Chebyshev  collocation  spectral  method  on  a  grid  of  9  points 
is  designed  to  differentiate  9  functions  exactly:  x’^  for  k  =  0, ...,  8.  Table  (3) 
is  designed  to  illustrate  the  same  two  points  at  Table  (2).  The  table  begins 
with  the  function  x®. 

Note  that  if  the  function  being  differentiated  is  not  exactly  an  integer 
frequency,  or  x^,  then,  say  for  Chebyshev,  differentiating  x®  ®  is  no  bet¬ 
ter  or  worse  than  the  result  for  differentiating  sin(5.5x).  The  point  that  is 
trying  to  be  made  is  that  the  differentiation  accuracy  of  spectral  collocation 
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Grid 

Pts 

Poly 

Order 

L2 

Error 

9 

2.2910-2® 

9 

3.5.5 

9.2610"^ 

9 

X® 

6.4910-2® 

9 

3.6.5 

to 

0 

1 

00 

9 

x" 

2.3910-2^ 

9 

1.7910-2 

9 

x8 

6.9610-24 

9 

x«-5 

4.2910-^ 

9 

x^ 

2.8310“ 

Table  3:  Chebyshev  Spectral  Collocation  Applied  to  Polynomials 


methods  on  a  finite  grid  of  size  N  is  accurate  with  the  accuracy  of  —  1, 
and  one  can  not  expect  that  a  wave-like  pulse  will  be  transmitted  better 
with  a  Fourier  spectral  method  than  with  a  Chebyshev  spectral.  The  only 
important  issue  is  the  dimension  of  the  space  and  the  boundary  conditions: 
use  Chebyshev  for  non-periodic  boundary  conditions  and  Fourier  for  period 
boundary  conditions. 

3.2  Very  High  Order  Finite  Differencing 

Now  suppose  that  build  a  series  of  algebraically  generated  finite  difference 
operators  of  increasing  accuracy  and  test  these  difference  operators  on  the 
function  sin(2x).  The  grid  size  will  be  fixed  at  33  points.  The  first  two  lines  of 
Table  (4)  illustrate  the  effect  of  the  Runge  phenomenon,  see  [5].  The  change 
in  the  error  from  periodic  boundary  conditions  on  a  uniform  grid  to  non¬ 
periodic  boundary  conditions  on  a  uniform  grid  is  from  10“^^  to  In 

addition,  note  that  applying  an  algebraic  polynomial  with  periodic  bound¬ 
ary  conditions  yields  a  result  comparable  to  applying  a  trigonometric,  i.e. 
Fourier  spectral  collocation,  polynomial  with  periodic  boundary  conditions. 
Furthermore,  one  can  observe  the  Runge  phenomenon  with  trigonometric 
polynomials  just  as  one  observes  it  with  algebraic  polynomials  when  the 
boundary  conditions  are  not  periodic.  Note  that  128  bit  arithmetic  is  being 
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■311 

Order 
of  Acc 

L2 

Error 

Periodic 

EC’s 

Grid  even 
or  Cheby 

32 

yes 

even 

— 

32 

no 

even 

33 

18 

7.4910"^^ 

no 

Cheby 

33 

20 

1.1710-1® 

64.0 

no 

Cheby 

33 

22 

1.7310-2° 

67.6 

no 

Cheby 

33 

24 

2.4310-22 

71.2 

no 

Cheby 

33 

26 

3.4910-2^ 

69.6 

no 

Cheby 

33 

28 

5.7510-2° 

60.7 

no 

Cheby 

33 

30 

1.3010-2^ 

44.2 

no 

Cheby 

33 

32 

8.8910-2° 

1.46 

no 

Cheby 

Table  4;  Finite  Difference  Accuracy  Approaching  Spectral  Accuracy 


used.  ^From  line  3  to  the  bottom  of  the  table  the  order  of  accuracy  is  in¬ 
creased  from  18  to  the  maximum,  i.e.,  spectral,  accuracy  of  32  is  obtained. 
When  one  tests  the  accuracy  of  a  finite  difference  operator  one  doubles  the 
grid  and  sees  the  error  decrease  as  (|)”  where  n  is  the  accuracy  of  the  scheme. 
This  comes  from  the  truncation  error  which  will  produce  a  factor  of  the  form 
(Ax)"^.  In  Table  (4),  it  is  the  number  n  which  is  being  increased  while  Ax 
remains  constant. 

Compare  Table  (4)  to  Table  (5)  in  which  a  Chebyshev  collocation  method 
on  an  increasing  grid  size  is  tested  on  sin(2x).  Note  that  in  the  following  table 
that  both  Ax  is  decreasing  and  n  is  increasing  in  the  expression  (Ax)”  as 
one  proceeds  down  the  table.  The  final  line  of  Table  (5)  is  the  Chebyshev 
method  on  a  grid  of  33  points  which  produces  a  result  compaxable  the  result 
in  Table  (4)  on  a  grid  of  33  points.  The  numbers  are  not  exactly  the  same 
because,  first  of  all,  all  calculations  are  near  machine  accuracy,  and,  second, 
the  differencing  coefficients  are  calculated  in  different  ways. 
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Grid 

Pts 

Alg 

Err 

Err 

Ratio 

9 

11 

13 

52.1 

15 

17 

2.3810-1° 

96.2 

19 

iSESiiii 

122.7 

21 

23 

3.0810-1° 

IfetelEM 

wm 

1.0910-2° 

33 

Table  5:  Chebyshev  Spectral  Collocation  of  Increasing  Order 
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3.3  Chebyshev  Spectral  Methods  and  Cosine  Polyno¬ 
mials 

This  subsection  will  review  Chebyshev  spectral  methods  and  the  equivalency 
with  cosine  polynomials.  Chebyshev  approximation  can  be  seen  as  approxi¬ 
mation  by  algebraic  polynomials, 

Toix)  =  1,  (42) 

Ti{x)  =  X, 

T2{x)  =  2x^  —  1, 

T3{x)  =  ix^  —  3x, 

or  as  approximation  by  a  cosine  series,  see  [7], 

n  n 

Tn{x)  =  cos(n  arccos  x)  =  cos{n0)  =  ^  ag(cos  Oy  =  ^  Qqx'^,  (43) 

g=0  9=0 

for  some  set  {a,}  and  where  x  =  cos{9).  K  one  now  chooses  a  numerical  grid 
defined  as  Xj  =  cos(j^),  j  =  0,  then  one  obtains  Tn{xj)  =  cos(^)  and 
the  pseudospectral  Chebyshev  method,  see  [10],  [21],  and  [2]. 

A  Chebyshev  spectral  method  involves  approximating  a  function,  /(a:), 
by  interpolating  an  algebraic  polynomial  to  point  values  f{xi)  where  the  grid 
points  are  given  by  the  ‘Chebyshev’  grid  points  x,  =  cos(0j).  An  equivalent 
process  is  to  interpolate  a  cosine  polynomial  to  the  evenly-spaced  point  values 
of  the  angle  di  and  to  consider  /(x)  to  be  evaluated  on  the  uniform  grid  of 
the  angle  variable  6i,  f(xi)  becomes  f{0i).  That  is,  see  [10],  if  the  Chebyshev 
series  for  /(x)  is 

OO 

Pf{x)  =  9{x)  =  (^kTk{x),  (44) 

A:=0 

then  the  expansion  coefficients  {a*}  can  be  found  in  two  equivalent  ways, 

ak  =  —  t  f{x)Tk{x){\  -  x'^y^'^dx  =  —  r  /(cos  0)  cos  k0d0  (45) 

By  this  transformation  of  the  independent  variable  one  can  perform  Cheby¬ 
shev  spectral  methods  on  a  uniform  grid  or  one  can  build  arbitrarily  high 
difference  operators  on  a  uniform  grid  which  have  stability  characteristics 
equivalent  to  the  usual  Chebyshev  spectral  method. 
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3.4  High- Order  Differencing  on  Chebyshev  Grids 

Chebyshev  spectral  methods  work  very  well  for  non-periodic  problems  pre¬ 
cisely  because  the  truncation  error  for  a  Chebyshev  polynomial  is  equil-ripple. 

As  shown  above,  the  truncation  error  for  polynomial  approximation, 
p„(x),  of  a  function  f{x)  is 

m  -  ?»(*)  =  (46) 

where  (  lies  between  the  smallest  and  the  largest  Xj.  If  one  wants  to  minimize 
the  error  due  to  the  term 


(x  —  Xo){x  -  Xi)...(x  —  Xn) 

then  the  n+1  sample  points  should  be  chosen  as  the  zeros  of  Chebyshev  poly¬ 
nomial  Tn+i{x),  see  [7].  This  selection  of  grid  points  ensures  that  the  error 
has  the  equal-ripple  property  which  is  characteristic  of  Chebyshev  polynomi¬ 
als.  A  Chebyshev  spectral  method  interpolates  the  highest  order  polynomial 
possible  onto  the  n  -b  1  degrees-of-freedom  defined  by  the  point  values  of  a 
function  at  the  n  -|-  1  zeros  of  Tn+i{x).  The  question  to  be  addressed  now 
is  what  if  the  grid  is  defined  as  the  zeros  of  Tn+x{x)  but  the  polynomial  is 
of  lower  order.  That  is,  n  could  be,  say,  128  whereas  the  polynomial  on 
this  grid  could  be  of  order  16.  It  is  well-known  that  high  order  polynomial 
interpolation  on  uniform  grids  is  essentially  an  ill-posed  problem. 

4  Wavelet-based  Grid  and  Order  Selection 

The  previous  section  introduced  the  idea  of  building  very  high  order  algebraically- 
generated  difference  operators  on  Chebyshev  grids  as  a  way  of  obtaining  very 
high  accuracy  which  is  almost  spectral  in  nature.  This  section  will  explore 
the  idea  of  performing  wavelet-based  grid  refinement  on  these  Chebyshev 
grids  as  a  way  to  obtain  the  necessary  Chebyshev  grid  distribution  near  a 
boundary  while  having  the  ability  to  refine  the  grid  away  from  the  boundary 
for  proper  physical-space  function  resolution. 
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4.1  A  Short  Review  of  Wavelets 

To  define  Daubechies-based  wavelets,  see  [6]  for  the  original  work  and  see 
[19]  for  an  introduction  to  wavelet-based  signal  processing,  consider  the  two 
functions  (j>(x),  the  scaling  function,  and  the  wavelet.  The  scaling 

function  is  the  solution  of  the  dilation  equation, 

L-l 

4>{x)  =  ^  hk<f>{2x  -  k),  (47) 

k=0 

where  ^{x)  is  normalized  <f>{x)dx  =  1,  and  the  wavelet  il){x)  is  defined  in 
terms  of  the  scaling  function, 

L-l 

V>(x)  =  y/2'^  gk^{2x  -  k).  (48) 

k=0 

One  builds  an  orthonormal  basis  from  4>{x)  and  ip{x)  by  dilating  and 
translating  to  get  the  following  functions: 

(f>{(x)  =  2~^(f>{2~^x  —  k),  (49) 

and 

V’i(x)  =  2-^{2-^x  -  k),  (50) 

where  j,  k  £  Z.  j  is  the  dilation  parameter  and  k  is  the  translation  pa¬ 
rameter.  The  coefficients  H  =  {Afclfcro  ^  —  {9k}kZl  a^re  related  by 
Qk  =  {—Vj^hL-k  for  k  =  0,...,!/  —  1.  All  wavelet  properties  are  specified 
through  the  parameters  H  and  G.  If  one’s  data  is  defined  on  a  continuous 
domain  such  as  f{x)  where  x  6  i?  is  a  real  number  then  one  uses  ^'[(x)  and 
ipki^)  fo  perform  the  wavelet  analysis.  If,  on  the  other  hand,  one’s  data  is 
defined  on  a  discrete  domain  such  as  f{i)  where  i  is  an  integer  then  the 
data  is  analyzed,  or  filtered,  with  the  coefficients  H  and  G.  In  either  case, 
the  scaling  function  ^(x)  and  its  defining  coefficients  H  detect  localized  low 
frequency  information,  i.e.,  they  are  low-pass  filters  (LPF),  and  the  wavelet 
^(x)  and  its  defining  coefficients  G  detect  localized  high  frequency  informa¬ 
tion,  i.e.,  they  are  high-pass  filters  (HPF).  Specifically,  H  and  G  are  chosen 
so  that  dilations  and  translations  of  the  wavelet,  form  an  orthonormal 

basis  of  L‘^{R)  and  so  that  ^(x)  has  M  vanishing  moments  which  determines 
the  accuracy.  In  other  words,  0fc(a:)  will  satisfy 

/OO 

(51) 

-OO 
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where  8ki  is  the  Kronecker  delta  function,  and  the  accuracy  is  specified  by 
requiring  that  =  V’o(®)  satisfy 

/CO 

'tp{x)x”^dx  =  0,  (52) 

-oo 

for  m  =  0, M  —  1.  Under  the  conditions  of  the  previous  two  equations,  for 
any  function  f{x)  €  L^{R)  there  exists  a  set  {djk}  such  that 

f{^)  =  a  djki/^Kx),  (53) 

j^Z  k^Z 

where  _ 

djk  =  /  f{x)ipl{x)dx.  (54) 

The  two  sets  of  coefficients  H  and  G  are  known  as  quadrature  mirror 
filters.  For  Daubechies  wavelets  the  number  of  coefficients  in  H  and  G,  or 
the  length  of  the  filters  H  and  G,  denoted  by  L,  is  related  to  the  number  of 
vanishing  moments  M  by  2M  =  L.  For  example,  the  famous  Haar  wavelet 
is  found  by  defining  H  as  ho  =  hi  =  1.  For  this  filter,  H,  the  solution  to 
the  dilation  equation  (47),  ^{x),  is  the  box  function:  4{x)  =  1  for  a;  G  [0,1] 
and  (f>{x)  =  0  otherwise.  The  Haar  function  is  very  useful  as  a  learning 
tool,  but  because  of  its  low  order  of  approximation  accuracy  and  lack  of 
differentiability  it  is  of  limited  use  as  a  basis  set.  The  coefficients  H  needed 
to  define  compactly  supported  wavelets  with  a  higher  degree  of  regularity 
can  be  found  in  [6] .  As  is  expected,  the  regularity  increases  with  the  support 
of  the  wavelet.  The  usual  notation  to  denote  a  Daubechies-based  wavelet 
defined  by  coefficients  H  of  length  L  is  Di. 

It  is  usual  to  let  the  spaces  spanned  by  (j>k{^)  parameter 

fe,  with  j  fixed,  be  denoted  by  Vj  and  Wj  respectively. 


Vj  =  4(x), 

(55) 

Wj  =  i^iix). 

^  kez  ^ 

(56) 

The  spaces  Vj  and  Wj 

are  related  by. 

...  C  Ui  C  Uo  C  V-i  c ..., 

(57) 
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and 


0^0+1,  (58) 

where  the  notation  Vo  =  Vi^Wi  indicates  that  the  vectors  in  Vi  are  orthog¬ 
onal  to  the  vectors  in  Wi  and  the  space  Vq  is  simply  decomposed  into  these 
two  component  subspaces. 

The  previously  stated  condition  that  the  wavelets  form  an  orthonormal 
basis  of  L'^{R)  can  now  be  written  as, 

L^R)  =  0  W,. 

jez 

Two  final  properties  of  the  spaces  Vj  are  that, 

n  =  {0}. 

i€Z 

and  _ 

U  =  L^R). 

j€Z 

4.2  Grid  Refinement  on  Uniform  Grids 

The  idea  of  using  wavelets  to  generate  numerical  grids  began  with  the  ob¬ 
servation  in  [12]  that  the  essence  of  an  adaptive  wavelet- Galerkin  method 
is  nothing  more  than  a  finite  difference  method  with  grid  refinement.  So, 
instead  of  letting  the  magnitude  of  wavelet  coefficients  choose  which  basis 
functions  to  use  in  a  Galerkin  approach,  let  the  same  coefficients  choose 
which  grid  points  to  use  and  then  think  of  the  wavelet  method  in  a  colloca¬ 
tion  sense. 

In  other  words,  suppose  a  calculation  begins  with  N  evenly-spaced  sam¬ 
ples  of  a  function  /  and  that  some  quadrature  method  produces  N  scaling 
function  coefficients  on  the  finest  scale  denoted  by  Vq.  If  the  spacing  between 
adjacent  values  in  the  vector  /  is  Ax  then  this  is  also  the  physical-space 
resolution  of  any  calculation  done  in  Vq.  Now,  decompose  Vq  once  to  get 
Vq  =  Vx®Wx.  Similarly  speaking,  the  physical  space  resolution  of  Vi  is 
2Ax  and  the  refinement  from  the  2Ax  physical-space  resolution  to  the  Ax 
physical-space  resolution  is  dictated  by  the  wavelet  coefficients  in  W\.  This 
is  the  reasoning  which  led  to  WOFD  and  to  the  following  subroutine  which  is 
at  the  heart  of  WOFD.  The  remainder  of  the  paper  is  concerned  with  giving 
the  reader  an  idea  of  how  the  WOFD  grid  refinement  software  works. 


(59) 

(60) 
(61) 
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4.3  Grid  Refinement  on  Chebyshev  Grids 

Chebysliev  grids  are  not  evenly-spaced  in  physical  space,  but  are  evenly- 
spaced  in  angle.  That  is,  a  Chebyshev  grid  comes  from  Xj  =  cos(0j),  j  = 
0,  where  the  angle  9j  =  ^  is  evenly-spaced.  The  above  described 

refinement  mechanism  can  now  be  applied  to  the  uniform  angle  grid  point 
values  to  define  a  new  numerical  grid.  That  is,  all  the  above  grid  refinement 
machinery  can  be  applied  to  Chebyshev  grids  where  each  subspace  Vj  will 
coincide  with  a  uniform  angle  or  usual  Chebyshev  grid  and  each  refinement 
subspace  Wj  will  coincide  with  additional  points  being  added  to  the  usual 
Chebyshev  grid.  It  is  well  known  that  the  Chebyshev  grid  is  the  best  grid,  in 
terms  of  minimal  error,  for  algebraic  polynomial  interpolation.  A  refinement, 
Wi,  on  Chebyshev  grid,  to  get  lo  =  ki  0  Wx  is  disigned  to  begin  with 
a  grid  which  in  one  sense  in  perfect,  the  Chebyshev  grid,  and  perturb  from 
this  grid. 


5  A  New  Numerical  Method:  WOFD2 

In  [15]  a  numerical  method  was  defined  which  was  called  the  Wavelet-Optimized 
Finite  Difference  method  or  WOFD.  WOFD  used  wavelets  in  their  finite 
difference  form.  In  essence,  this  meant  that  wavelets  were  used  to  choose 
a  numerical  grid  and  all  computations  were  performed  on  this  grid  using 
arbitrary-grid  finite  difference  operators. 

WOFD2  is  an  extension  of  this  idea.  WOFD2  uses  wavelets  to  choose  not 
only  a  numerical  grid  but  also  the  order  of  the  difference  operator  used  on  this 
grid.  In  addition,  WOFD2  uses  very  high  order  finite  difference  operators 
on  the  order  of  8,  16  or  maybe  32.  Furthermore,  the  physical-space  grids 
are  no  longer  evenly-spaced  at  every  resolution  but  are  Chebyshev.  That  is, 
wavelet-based  grid  generation,  see  [16],  requires  that  a  grid  be  selected  from 
a  uniform  finest  grid.  But,  high  order  polynomials  can  be  highly-oscillatory 
on  uniform  grids.  Therefore,  WOFD2  works  with  Chebyshev  grids  at  each 
resolution  level.  Recall,  that  Chebyshev  grids  Xi  =  cos(^,)  are  not  uniform 
in  the  physical  space  variable  x,  but  are  uniform  in  the  angle  variable  6i.  It 
is  in  this  uniform  angle  variable  6i  that  grid  refinement  is  performed. 

Using  the  grid  selection  mechanism  in  [16]  to  select  order  is  a  minor  exten¬ 
sion  of  the  idea  that  wavelets  are  very  good  at  finding  regions  of  the  domain 
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at  which  a  large  numerical  error  is  likely  to  occur.  Numerical  error  will  deter¬ 
mined  by  the  truncation  error  of  a  polynomial  which  is  locally  interpolated 
to  the  data.  The  truncation  error  will  be  the  product  of  intervals  and  a  con¬ 
stant.  Imagine  the  intervals  are  all  a  multiple  of  a  smallest  interval  Ax  then 
the  key  component  in  the  truncation  error  will  be  (Ax)”.  This  component 
can  be  decreased  by  either  decreasing  the  size  of  Ax  or  by  increasing  the 
order  of  the  scheme,  i.e.,  increasing  n.  Or,  one  can  decrease  Ax  and  increase 
n  simultaneously. 

5.1  Comparison  with  hp-Refinement 

In  the  finite  element  literature,  see  [20],  [1],  [11],  the  idea  of  refining  the  grid 
and  increasing  the  polynomial  order  is  known  as  hp-refinement.  The  theory 
from  hp-refinement  can  certainly  be  applied  to  the  new  method  WOFD2, 
even  though  WOFD2  works  with  polynomials  only  for  generation  of  finite 
difference  operators  to  be  applied  in  the  physical  space.  One  of  the  most 
important  results  from  hp-refinement  theory  from  which  WOFD2  can  ben¬ 
efit  is  that  when  the  function  being  differentiated  is  smooth  then  the  rate 
of  convergence  is  controlled  by  the  polynomial  degree.  For  the  purpose  of 
pulse  propagation  in  aeroacoustics  it  is  apparent  that  a  high  order  differen¬ 
tiation,  i.e.,  high  order  polynomial  interpolation,  will  propagate  the  pulse 
more  faithfully  than  grid  refinement  on  the  same  pulse,  assuming  the  pulse 
is  smooth. 

5.2  Numerical  Experiments  with  WOFD2 

This  section  will  provide  the  results  from  numerous  numerical  experiments 
performed  with  WOFD2.  For  all  the  numerical  experiments  in  this  section  of 
the  paper  a  Gaussian  pulse  enters  the  domain  from  the  right-hand  side  and 
travels  to  the  left.  The  governing  equation  is  the  1  dimensional  hyperbolic 
wave  equation, 

Ut{x,t)  =  U,{x,t),  i7(x,0)=e-‘^(*-)^  (62) 

for  some  constant  c.  See  Figure  1  for  a  plot  of  this  initial  condition.  Note  that 
the  domain  extends  from  0  to  jt.  The  final  time  for  all  simulation  is  7r/2.  The 
simulation  is  stopped  at  this  value  because  at  7r/2  the  Chebyshev  grid  has  a 
maximum  spacing  between  grid  points  and  hence  a  minimum  resolution. 
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Grid 

Pts 

Order 
of  Acc 

i/2 

Error 

L<x> 

Error 

Final 

Time 

64 

8 

ismiiBn 

7r/2 

64 

16 

BglQESi 

7r/2 

64 

32 

7.10  ♦  10-® 

7r/2 

64 

48 

7r/2 

Table  6:  W0FD2  of  Accuracies  8,  16,  32,  and  48 


5.2.1  No  Adaptation 

First  we  consider  the  case  of  very  high  order  finite  differencing  on  a  Cheby- 
shev  grid.  The  grid  size  is  kept  fixed  at  64  points,  and  the  order  is  increased 
from  8  to  48.  The  errors  decrease  in  a  nice  and  uniform  manner.  No  unusual 
numerical  oscillations  occur. 


5.2.2  Adapting  Grid  Only;  Order  8  Spatial  Differencing 

In  this  subsection  the  order  of  the  spatial  differencing  is  kept  fixed  at  order  8. 
No  results  are  found  for  threshold  values  of  10“^  and  10“^.  This  is  because  it 
seems  to  be  a  characteristic  of  adaptive  methods  that  a  very  rough  threshold 
value  can  degrade  the  performance  of  the  method.  It  is  better  to  start  with 
threshold  values  less  than  or  equal  to  10“^.  The  first  row  of  the  table  is  the 
worst  possible  performance  where  no  refinement  is  done  and  grid  is  64  points, 
and  the  last  row  of  the  table  is  best  possible  performance  where  the  grid  is  128 
points.  Note  that  the  software  is  constructed  to  work  with  both  periodic  and 
non-periodic  boundary  conditions,  so  that  when  the  boundary  conditions  are 
non-periodic  the  possible  number  of  points  becomes  2^-f  1  which  includes  the 
right-hand  boundary  point.  For  periodic  boundary  conditions  the  number  of 
grid  points  is  2^  since  the  right-hand  boundary  point  is  equal  to  the  first 
point  on  the  left-hand  boundary. 
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Grid 

Pts 

Grid 

Order 
of  Acc 

L2 

Error 

Loo 

Error 

Thresh 

Final 

Time 

128/64 

65 

8 

3.82  =«  10"'' 

1.78*10-2 

100.0 

7r/2 

128/64 

77 

8 

1.67*10"^ 

8.66  *  lO-"" 

10-3 

7r/2 

128/64 

79 

8 

1.63  *  10-^ 

8.76  *  10"^ 

10"^ 

7r/2 

128/64 

82 

8 

8.24  *  10-® 

3.24  *  10-^ 

10-3 

7r/2 

128/64 

84 

8 

8.07  *  10-® 

3.24  *  10-^ 

10-3 

7r/2 

128 

129 

8 

6.51  *  10-® 

3.24  *  10-^ 

0.0 

7r/2 

Table  7:  W0FD2  Adaptive  Grid,  but  Accuracy  fixed  at  8 


5.2.3  Adapting  Grid  Only:  Order  16  Spatial  Differencing 

Much  of  what  was  said  for  the  8th  order  table  above  can  be  said  here.  The 
second  row  of  the  following  table  shows  how  the  performance  can  be  slightly 
degraded  for  relatively  large  threshold  values.  In  this  case  the  degradation 
occurs  at  the  threshold  value  of  10“^.  This  is  a  minor  point.  Generally 
speaking,  just  start  with  a  smaller  threshold  value. 

5.2.4  Adapting  the  Grid  and  Order 

This  final  table  is  the  culmination  of  the  paper  and  an  example  of  WOFD2 
with  all  options  in  use.  The  grid  is  adjusted  between  a  maximum  density  of 
128  and  a  minimum  density  of  64.  The  order  of  accuracy  is  adjusted  between 
a  maximum  order  of  16  and  a  minimum  order  of  8.  The  error  converge  in  a 
nice  manner  toward  the  minimum  error  which  occurs  at  the  maximum  grid 
density  of  128  and  the  maximum  order  of  accuracy  of  16. 

The  usual  Chebyshev  grid  is  evenly-spaced  in  angle  6i  =  iir/N  for  i  = 
0,...,N.  In  the  physical  space  the  grid  distribution  is  Xi  =  cos(0,)  which 
is  shaped  like  a  semi-circle.  When  one  applies  the  wavelet  grid  adaptation 
to  this  evenly-spaced  $i  then  obtains  in  the  physical  space  the  distribution 
Xi  =  cos{6i)  in  the  portion  of  the  domain  away  from  the  pulse,  and  the  twice- 
as-dense  grid  distribution  Xi  =  cos^iir / {2N))  in  the  portion  of  the  domain 
near  the  pulse.  Note  that  the  grid  is  the  usual  Chebyshev  grid  near  the 
boundary.  It  is  only  safely  away  from  the  boundary  that  the  grid  density 
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Grid 

Pts 

Grid 

Order 
of  Acc 

1/2 

Error 

Loo 

Error 

Thresh 

Final 

Time 

128/64 

65 

16 

ISjiQSjgil 

100.0 

x/2 

128/64 

81 

16 

■B 

128/64 

16 

■B 

^@1 

89 

16 

8.09  ♦  10-5 

10-5 

7r/2 

87 

16 

9.20  + 10-5 

10-5 

x/2 

128/64 

91 

16 

BOB! 

wnM 

128/64 

93 

16 

128 

129 

16 

4.26  *  10-« 

0.0 

1^12 

Table  8:  W0FD2  Adaptive  Grid,  but  Accuracy  fixed  at  16 


Grid 

Density 

Grid 

Order 
of  Acc 

Z/2 

Error 

Loo 

Error 

Thresh 

Final 

Time 

64 

65 

8 

100.0 

x/2 

128/64 

81 

HSi 

x/2 

10-^ 

■B 

IB&ESi 

mm 

■B 

84 

16/8 

10-5 

x/2 

86 

16/8 

6.49  =<=  10-"' 

10-^ 

x/2 

128/64 

87 

16/8 

BBKiii 

BBUsBI 

10-5 

128/64 

90 

16/8 

5.12=1=10-5 

BSHKiBI 

10-^ 

128 

129 

16 

ibiiqs 

0.0 

x/2 

Table  9:  W0FD2  with  Grid  and  Order  Adaptation 
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makes  an  abrupt  change  in  density.  See  Figure  2  for  an  example  of  an  initial 
grid. 

If  the  numerical  scheme  is  working  properly  then  the  pulse  will  propagate 
to  the  middle  of  the  domain  and  be  similar  in  shape  to  the  initial  condition. 
The  best  measure  of  this  similarity  is  the  Too  error.  At  the  final  time  the 
pulse  will  appear  as  in  Figure  3. 

The  Chebyshev  grid  is  naturally  more  dense  near  the  boundaries  than  in 
the  middle  of  the  domain.  With  the  wavelet  adaptation  of  this  Chebyshev 
grid,  the  grid  points  can  be  kept  dense  while  maintaining  a  the  Chebyshev 
distribution  throughout  most  of  the  domain.  Again,  the  most  important 
region  of  the  domain  for  a  Chebyshev  distribution  is  near  the  boundary.  See 
Figure  4  for  the  grid  distribution  at  the  final  time  when  the  pulse  has  reached 
the  middle  of  the  domain. 

Without  grid  refinement  or  order  refinement  the  peak  numerical  error  at 
the  final  time  should  be  near  the  peak  of  the  pulse,  since  it  is  this  portion  of 
the  function  which  is  most  difficult  to  represent  by  polynomial  interpolation. 
See  Figure  (5)  for  an  example  of  such  an  error. 

If  the  wavelet  refinement  threshold  is  not  sufficiently  low  then  one  will 
see  the  peak  error  appear  near  a  region  of  the  domain  where  there  is  a  grid 
or  stencil  discontinuity.  A  ‘sufficiently  low’  refinement  threshold  will  on  the 
order  of  the  L<x>  error  when  no  grid  or  order  refinement  is  executed.  In 
Figure  6  noise  is  amplified  at  the  interface  where  both  the  stencil  and  grid 
are  refined.  If  the  wavelet  refinement  threshold  is  adjusted  to  a  smaller  value 
then  one  can  obtain  an  error  profile  similar  to  that  in  Figure  (5). 

When  both  the  stencil  and  grid  are  changed  throughout  the  calculation, 
one  finds  a  relatively  wide  stencil  near  the  peak  value  of  the  pulse.  For  the 
example  of  a  17  point  stencil  with  accuracy  of  16  at  the  pulse  and  a  9  point, 
accuracy  8,  stencil  away  from  the  pulse  see  Figure  (7). 

6  Conclusion 

This  paper  has  covered  many  topics  related  to  the  construction  of  a  very 
high  order  adaptive  order  and  adaptive  grid  numerical  method  which  has 
been  named  the  Wavelet-Optimized  Difference  Method  2,  or  W0FD2.  First 
it  was  necessary  to  explore  the  various  ways  in  which  difference  operators  can 
be  constructed.  This  included  a  comparison  of  difference  operators  generated 
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from  algebraic,  trigonometric,  exponential,  and  cosine  polynomials.  Next, 
wMch  type  of  polynomial  would  be  best  for  tbe  construction  of  very  higb 
order  numerical  differencing.  Tbe  conclusion,  which  is  not  a  big  surprise, 
is  that  one  should  use  algebraic  polynomials  on  Chebyshev  grids.  The  next 
step  was  to  apply  wavelet  grid  and  order  adaptation  in  order  to  be  able 
to  reduce  errors  throughout  the  domain  by  either  increasing  the  order  of 
the  numerical  method  or  by  increasing  the  grid  density  in  the  appropriate 
region.  The  results  of  the  numerical  tests  were  very  positive  and  it  appears 
that  WOFD2  will  applicable  to  a  large  range  of  numerical  problems.  The 
version  of  WOFD2  which  has  been  presented  here  has  been  ‘tweaked’  very 
little.  That  is,  it  worked  essentially  for  the  first  time  it  was  tried.  This  is 
encouraging  because  most  high  order  numerical  methods  require  some  kind  of 
filtering  or  other  refinement.  Future  plans  for  WOFD2  would  be,  perhaps,  to 
try  to  find  a  proof  of  stability  and  to  apply  the  method  in  higher  dimensions. 


29 


References 

[1]  I.  Babuska,  “The  p-  and  hp- Versions  of  the  Finite  Element  Method: 
the  State  of  the  Art.”  Finite  Elements:  Theory  and  Applications^  edited 
by  D.L.  Dwoyer,  M.Y.  Hussaini  and  R.G.  Voigt.  New  York:  Springer- 
Verlag.  pp.  199-239  (1988). 

[2]  C.  Canuto,  M.Y.  Hussaini,  A.  Quarteroni,  T.A.  Zang,  (1988)  “Spectral 
Methods  in  Fluid  Dynamics” ,  Springer- Verlag. 

[3]  M.  Carpenter  and  D.  Gottlieb,  “Spectral  Methods  on  Arbitrary  Grids”, 
ICASE  Report  No.  95-37,  NASA  CR-198158. 

[4]  S.N.  Christofi,  “The  Study  of  Building  Blocks  for  Essentually  Non- 
Oscillatory  (ENO)  Schemes”,  PhD  Thesis,  Division  of  Applied  Math¬ 
ematics,  Brown  University,  May  1996. 

[5]  G.  Dahlquist  and  A.  Bjorck,  (1974)  “Numerical  Methods”,  Prentice- 
Hall. 

[6]  I.  Daubechies,  (1988)  “Orthonormal  Basis  of  Compactly  Supported 
Wavelets”,  Comm.  Pure  Appl.  Math.,  41  pp.  909-996. 

[7]  P.J.  Davis,  (1975)  “Interpolation  and  Approximation”,  Dover. 

[8]  W.S.  Don  and  C.B.  Quillen,  “Numerical  Simulation  of  Shock- Cylinder 
Interactions”,  Journal  of  Computational  Physics,  122,  244-265  (1995). 

[9]  (1996)  G.  Erlebacher,  M.Y.  Hussaini,  L.  Jameson,  “Wavelets:  Theory 
and  Applications”,  Oxford. 

[10]  D.  Gottlieb  and  S.A.  Orszag,  (1977)  “Numerical  Analysis  of  Spectral 
Methods:  Theory  and  Applications”,  SIAM-CBMS,  Philadelphia. 

[11]  W.  Gui  and  I.  Babuska,  (1986),  “The  h-,  p-  and  hp- Versions  of  the 
Finite  Element  Method  in  One  Dimension.  Part  I:  The  Error  Analysis 
of  the  p- Version.  Part  II:  The  Error  Analysis  of  the  h  and  hp- Versions. 
Part  III:  The  Adaptive  hp- Version.”  Numerische  Mathematik,  Vol.  49, 
pp.  577-612;  613-657;  659-683. 


30 


[12]  L.  Jameson,  “On  The  Wavelet  Based  Differentiation  Matrix”,  Journal 
of  Scientific  Computing,  September  1993  and  ICASE  Report  No.  93-95, 
NASA  CR-191583. 

[13]  L.  Jameson,  “On  The  Differentiation  Matrix  for  Daubechies-Based 
Wavelets  on  an  Interval”,  SIAM  J.  Sci.  Comp.,  Vol.  17,  Issue  2,  3/96 
and  ICASE  Report  No.  93-94,  NASA  CR-191582. 

[14]  L.  Jameson,  T.L.  Jackson,  D.G.  Lcisseigne,  “Wavelets  as  a  Numerical 
Tool”,  Proceedings  from  the  Joint  US- Japan  Workshop  on  Combustion, 
1993,  J.  Buckmaster,  T.  Takeno  (eds.).  Springer- Verlag. 

[15]  L.  Jameson,  “On  the  Wavelet- Optimized  Finite  Difference  Method”, 
ICASE  Report  No.  94-9,  NASA  CR-191601. 

[16]  L.  Jameson,  “Wavelet-based  Grid  Generation”,  12-95,  submitted  to 
Journal  of  Computational  Physics. 

[17]  G.E.  Karniadakis  and  S.A.  Orszag,  Physics  Today  34,  (1993). 

[18]  H.O.  Kreiss,  technical  report.  University  of  Uppsala,  Sweden,  1978. 

[19]  G.  Strang  and  T.  Nguyen,  “Wavelets  and  Filter  Banks”,  Wellesley- 
Cambridge  Press,  1996. 

[20]  B.  Szabo  and  I.  Babuska,  (1991)  “Finite  Element  Analysis”,  Wiley  In¬ 
terscience. 

[21]  R.G.  Voigt,  D.  Gottlieb,  Y.  Hussaini,  (1984)  “Spectral  Methods  for  Par¬ 
tial  Differential  Equations”,  SIAM. 


31 


32 


Figure  2:  Initial  Grid  Density  For  Pulse  Entering  Domain 
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